Universaltool Regressionsanalyse II (Mehrebenenmodelle)

Tag 3: Using Multi-Level Models

Samuel Merk

08.04.2022

Materials von Tag 1 und 2

Die Materials von Tag 1 & 2 könnt ihr hier herunterladen. Alle files sind auf unter https://bit.ly/merk039, die Dokumentation des heutigen Tag 3 findet ihr unter https://bit.ly/merk035.

Mein Plan für Tag III

  • Wdh: Einfache Multi-Level Begriffe konzeptuell verstehen
    • No Pooling, Complete Pooling, Shrinkage
    • Intercept Only Model, Random Intercept Model, Random Slope Model, Random Intercept and Random Slope Model
    • Intraklassenkorrelation
    • Fixed Effects, Random Effect
  • Unsere Beispieldaten
  • Einfache Multi-Level Modelle parametrisieren und schätzen

Unsere Beispieldaten

Sleep Study

Im Paket {lme4} ist der Datensatz sleepstudy enthalten. Diesen Daten liegt ein Experiment zugrunde indem dem Pbd Schlaf entzogen wurde Für uns wichtige Variablen sind:

  • Reaction: Durchschnittliche Reaktionszeit in einem Reaktionszeitexperiment in (ms)
  • Days: Anzahl der Tage des Schlafentzugs (der Schlafentzug startete nach dem zweiten Tag)
  • Subject: Personidentifier
library(lme4)
head(sleepstudy)
Reaction Days Subject
249.5600 0 308
258.7047 1 308
250.8006 2 308
321.4398 3 308
356.8519 4 308
414.6901 5 308

BilaKi Daten

data_bilaki_wide <- 
  read_sav("data/BilaKi_wide.sav") 

head(data_bilaki_wide)
ID V_ID Tgroup m0_s1 m1_s1 m2_s1
1 1 1 2.166667 3.166667 NA
2 1 1 2.333333 NA 3.166667
3 1 1 3.000000 3.416667 NA
4 1 1 3.666667 NA NA
5 1 1 3.250000 4.000000 4.500000
6 1 1 4.500000 3.416667 4.833333

Elke hat uns Daten aus ihrem Projekt mitgebracht. Diese liegen im wide-Format vor und müssen zunächst einmal in das long-Format gebracht werden. Das kann bequem mit dem im tidyverse enthaltenen Paket {tidyr} geschehen:

data_bilaki_long <- 
  data_bilaki_wide %>% 
  pivot_longer(names_to = "variable", 
               values_to = "perf",
               cols = c(m0_s1, m1_s1, m2_s1))

head(data_bilaki_long)
ID V_ID Tgroup variable perf
1 1 1 m0_s1 2.166667
1 1 1 m1_s1 3.166667
1 1 1 m2_s1 NA
2 1 1 m0_s1 2.333333
2 1 1 m1_s1 NA
2 1 1 m2_s1 3.166667

Dann muss man noch aus den Variablennamen die Info über den Messzeitpunkt entnehmen und eine nominale Zeitvariable erstellen.

data_bilaki_long <- 
  data_bilaki_long %>% 
  mutate(time = substr(variable, 2, 2),
         time_factor = factor(time))

head(data_bilaki_long)
ID V_ID Tgroup variable perf time time_factor
1 1 1 m0_s1 2.166667 0 0
1 1 1 m1_s1 3.166667 1 1
1 1 1 m2_s1 NA 2 2
2 1 1 m0_s1 2.333333 0 0
2 1 1 m1_s1 NA 1 1
2 1 1 m2_s1 3.166667 2 2

Popularity Data (Hox et al., 2017)

Dieser Datensatz entstammt den einführenden Beispielen des gut lesbaren Lehrbuchs von Joop Hox et al (2017) Für uns wichtige Variablen sind:

  • popular: Eine Likertskala zur Selbsteinschätzung der Beliebtheit
  • extrav: Extraversion (Big Five)
  • texp: Berufserfahrung der Lehrkraft

Import

data_popularity <- read_sav(file ="https://github.com/MultiLevelAnalysis/Datasets-third-edition-Multilevel-book/blob/master/chapter%202/popularity/SPSS/popular2.sav?raw=true")

head(data_popularity)
pupil class extrav sex texp popular popteach Zextrav Zsex Ztexp Zpopular Zpopteach Cextrav Ctexp Csex
1 1 5 1 24 6.3 6 -0.1703149 0.9888125 1.486153 0.8850133 0.6690561 -0.215 9.737 0.5
2 1 7 0 24 4.9 5 1.4140098 -1.0108084 1.486153 -0.1276291 -0.0430845 1.785 9.737 -0.5
3 1 4 1 24 5.3 6 -0.9624772 0.9888125 1.486153 0.1616973 0.6690561 -1.215 9.737 0.5
4 1 3 1 24 4.7 5 -1.7546396 0.9888125 1.486153 -0.2722923 -0.0430845 -2.215 9.737 0.5
5 1 5 1 24 6.0 6 -0.1703149 0.9888125 1.486153 0.6680185 0.6690561 -0.215 9.737 0.5
6 1 4 0 24 4.7 5 -0.9624772 -1.0108084 1.486153 -0.2722923 -0.0430845 -1.215 9.737 -0.5

Wdh: Begriffe

AA: Elaboriert diese Begriffe in dem Ihr Sie auf eines der händisch skizzierten Datenbeispiele (siehe live-demo) anwendet.

Einfache Multi-Level Modelle parametrisieren und schätzen

Sleep Study

Die Mehrebenenstruktur verstehen

Fragestellung

Wie kann die Mehrebenenstruktur der Sleep Study Datensatzes beschrieben werden?

Lösungshinweis 1

Mit head(sleepstudy) erkennt man, dass der Datensatz nur aus AV, UV, Subjectidentifier besteht. Also muss können maximal zwei Ebenen im Datensatz vorkommen.
Bildet man ein Tabelle mit Days und Subject kann man entscheiden ob eine perfect iode rimperfecte Hierarchie vorliegt.

Lösungshinweis 2

Eine solche Tabelle erhält man mit table(sleepstudy$Days, sleepstudy$Subject)

Lösung

Es handelt um eine perfekte Zweiebenenstruktur mit folgenden Nestgrößen und folgender Nestanzahl:

library(lme4) #enthält den Datensatz
# Nestanzahl
length(unique(sleepstudy$Subject))
## [1] 18
# Nestgrößen
sleepstudy %>% 
  group_by(Subject) %>% 
  summarize(Nestgröße = n())
Subject Nestgröße
308 10
309 10
310 10
330 10
331 10
332 10
333 10
334 10
335 10
337 10
349 10
350 10
351 10
352 10
369 10
370 10
371 10
372 10

Intraklasenkorrelation schätzen und verstehen

Fragestellung

Der ICC hilft dabei empirisch zu erkennen wie stark sich die Mehrebenenstruktur in der abhängigen Variable ausdrückt. Dazu schätzt man ein Intercept-Only Modell und setzt die Intercept-Varianz ins Verhältnis zur Gesamtvarianz.

\[ \begin{aligned} \operatorname{Reaction}_{i} &\sim N \left(\alpha_{j[i]}, \sigma^2 \right) \\ \alpha_{j} &\sim N \left(\mu_{\alpha_{j}}, \sigma^2_{\alpha_{j}} \right) \text{, for Subject j = 1,} \dots \text{,J} \end{aligned} \]

Lösungshinweis 1

Die Syntax für das Intercept-Only Modell lautet lmer(Reaction ~ 1 + (1|Subject), data = sleepstudy)

Lösungshinweis 2

Richtig schnell geht es mit der Funktion icc() aus dem Package {performance}.

Lösung

library(performance)
library(lme4)
mod01_sleepstudy <- lmer(Reaction ~ 1 + (1|Subject), data = sleepstudy)
icc(mod01_sleepstudy)
## # Intraclass Correlation Coefficient
## 
##      Adjusted ICC: 0.395
##   Conditional ICC: 0.395

Dieser ICC ist vglw. groß! Multi-Level Modelling ist also nicht nur theoretisch impliziert sonder auch empirisch.

Single Level Model schätzen

Fragestellung

Um später entscheiden zu können inwiefern Random Interept und Slope empirisch gerechtfertigt sind, kann man zunächst ein Single Level Model schätzen um dieses später mit den komplexeren Multi-Level Modellen zu schätzen.

Lösung

mod02_sleepstudy <- lm(Reaction ~ Days, data = sleepstudy)
summary(mod02_sleepstudy)
## 
## Call:
## lm(formula = Reaction ~ Days, data = sleepstudy)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -110.848  -27.483    1.546   26.142  139.953 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  251.405      6.610  38.033  < 2e-16 ***
## Days          10.467      1.238   8.454 9.89e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 47.71 on 178 degrees of freedom
## Multiple R-squared:  0.2865, Adjusted R-squared:  0.2825 
## F-statistic: 71.46 on 1 and 178 DF,  p-value: 9.894e-15

Das entspricht einem Modell \[ \operatorname{Reaction} = {\color{#8cd000}{\beta_{0}}} + {\color{#8cd000}{\beta}}_{1}(\operatorname{Days}) + {\color{#8cd000}{\epsilon}} \] mit den Koeffizienten \[ \operatorname{\widehat{Reaction}} = 251.41 + 10.47(\operatorname{Days}) \]

Random Intercept Model schätzen und vergleichen

Fragestellung

Wie ändern sich die Koeffizienten des Modells, wenn man zulässt, dass Intercepts zwischen den Clustern variieren?

Lösungshinweis

Das Intercept kann ja als konstanter Prädiktor mit dem Wert 1 aufgefasst werden. Daher lautet die Syntax für das Random Intercept (1|Subject) (lies: »das Intercept darf über die Clustervariable Subject variieren«).

Lösung

mod03_sleepstudy <- lmer(Reaction ~ Days + (1|Subject), data = sleepstudy)
summary(mod03_sleepstudy)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Reaction ~ Days + (1 | Subject)
##    Data: sleepstudy
## 
## REML criterion at convergence: 1786.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2257 -0.5529  0.0109  0.5188  4.2506 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Subject  (Intercept) 1378.2   37.12   
##  Residual              960.5   30.99   
## Number of obs: 180, groups:  Subject, 18
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 251.4051     9.7467   25.79
## Days         10.4673     0.8042   13.02
## 
## Correlation of Fixed Effects:
##      (Intr)
## Days -0.371

Das entspricht einem Modell

## Warning: Colorization of greek notation not currently implemented for merMod
## models

\[ \begin{aligned} \operatorname{Reaction}_{i} &\sim N \left(\alpha_{j[i]} + \beta_{1}(\operatorname{Days}), \sigma^2 \right) \\ \alpha_{j} &\sim N \left(\mu_{\alpha_{j}}, \sigma^2_{\alpha_{j}} \right) \text{, for Subject j = 1,} \dots \text{,J} \end{aligned} \] mit den Koeffizienten

## Warning: Colorization of greek notation not currently implemented for merMod
## models

\[ \begin{aligned} \operatorname{\widehat{Reaction}}_{i} &\sim N \left(251.41_{\alpha_{j[i]}} + 10.47_{\beta_{1}}(\operatorname{Days}), \sigma^2 \right) \\ \alpha_{j} &\sim N \left(0, 37.12 \right) \text{, for Subject j = 1,} \dots \text{,J} \end{aligned} \] Die Modelle vergleichen kann man mit der Syntax

tab_model(mod01_sleepstudy, mod02_sleepstudy, mod03_sleepstudy)
  Reaction Reaction Reaction
Predictors Estimates CI p Estimates CI p Estimates CI p
(Intercept) 298.51 280.65 – 316.37 <0.001 251.41 238.36 – 264.45 <0.001 251.41 232.17 – 270.64 <0.001
Days 10.47 8.02 – 12.91 <0.001 10.47 8.88 – 12.05 <0.001
Random Effects
σ2 1958.87   960.46
τ00 1278.34 Subject   1378.18 Subject
ICC 0.39   0.59
N 18 Subject   18 Subject
Observations 180 180 180
Marginal R2 / Conditional R2 0.000 / 0.395 0.286 / 0.282 0.280 / 0.704

Die Modelle gegeneinander Testen \(H_0: \text{beide Modelle klären gleich viel Varianz auf}\) kann man mit dem anova()-Befehl:

anova(mod03_sleepstudy, mod02_sleepstudy)
## refitting model(s) with ML (instead of REML)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
mod02_sleepstudy 3 1906.293 1915.872 -950.1465 1900.293 NA NA NA
mod03_sleepstudy 4 1802.079 1814.850 -897.0393 1794.079 106.2144 1 0

Unterschied ist signifikant, also steigert das Random Intercept die Modellpassung überzufällig.

Random Intercept Random Slope Model schätzen und vergleichen

Fragestellung

Wie ändern sich die Koeffizienten des Modells, wenn man zulässt, dass Intercept und Slopes variieren?

Lösungshinweis

Die Syntax für ddie Random Effects lautet (1 + Dayes|Subject) (lies: »das Intercept und die Slopes dürfen über die Clustervariable Subject variieren«).

Lösung

mod04_sleepstudy <- lmer(Reaction ~ Days + (1 + Days|Subject), data = sleepstudy)
summary(mod03_sleepstudy)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Reaction ~ Days + (1 | Subject)
##    Data: sleepstudy
## 
## REML criterion at convergence: 1786.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2257 -0.5529  0.0109  0.5188  4.2506 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Subject  (Intercept) 1378.2   37.12   
##  Residual              960.5   30.99   
## Number of obs: 180, groups:  Subject, 18
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 251.4051     9.7467   25.79
## Days         10.4673     0.8042   13.02
## 
## Correlation of Fixed Effects:
##      (Intr)
## Days -0.371

Das entspricht einem Modell

## Warning: Colorization of greek notation not currently implemented for merMod
## models

\[ \begin{aligned} \operatorname{Reaction}_{i} &\sim N \left(\alpha_{j[i]} + \beta_{1j[i]}(\operatorname{Days}), \sigma^2 \right) \\ \left( \begin{array}{c} \begin{aligned} &\alpha_{j} \\ &\beta_{1j} \end{aligned} \end{array} \right) &\sim N \left( \left( \begin{array}{c} \begin{aligned} &\mu_{\alpha_{j}} \\ &\mu_{\beta_{1j}} \end{aligned} \end{array} \right) , \left( \begin{array}{cc} \sigma^2_{\alpha_{j}} & \rho_{\alpha_{j}\beta_{1j}} \\ \rho_{\beta_{1j}\alpha_{j}} & \sigma^2_{\beta_{1j}} \end{array} \right) \right) \text{, for Subject j = 1,} \dots \text{,J} \end{aligned} \] mit den Koeffizienten

## Warning: Colorization of greek notation not currently implemented for merMod
## models

\[ \begin{aligned} \operatorname{\widehat{Reaction}}_{i} &\sim N \left(251.41_{\alpha_{j[i]}} + 10.47_{\beta_{1j[i]}}(\operatorname{Days}), \sigma^2 \right) \\ \left( \begin{array}{c} \begin{aligned} &\alpha_{j} \\ &\beta_{1j} \end{aligned} \end{array} \right) &\sim N \left( \left( \begin{array}{c} \begin{aligned} &0 \\ &0 \end{aligned} \end{array} \right) , \left( \begin{array}{cc} 24.74 & 0.07 \\ 0.07 & 5.92 \end{array} \right) \right) \text{, for Subject j = 1,} \dots \text{,J} \end{aligned} \] Die Modelle vergleicht man mit der Syntax

tab_model(mod01_sleepstudy, mod02_sleepstudy, mod03_sleepstudy, mod04_sleepstudy)
  Reaction Reaction Reaction Reaction
Predictors Estimates CI p Estimates CI p Estimates CI p Estimates CI p
(Intercept) 298.51 280.65 – 316.37 <0.001 251.41 238.36 – 264.45 <0.001 251.41 232.17 – 270.64 <0.001 251.41 237.94 – 264.87 <0.001
Days 10.47 8.02 – 12.91 <0.001 10.47 8.88 – 12.05 <0.001 10.47 7.42 – 13.52 <0.001
Random Effects
σ2 1958.87   960.46 654.94
τ00 1278.34 Subject   1378.18 Subject 612.10 Subject
τ11       35.07 Subject.Days
ρ01       0.07 Subject
ICC 0.39   0.59 0.72
N 18 Subject   18 Subject 18 Subject
Observations 180 180 180 180
Marginal R2 / Conditional R2 0.000 / 0.395 0.286 / 0.282 0.280 / 0.704 0.279 / 0.799

Die Modelle gegeneinander Testen \(H_0: \text{beide Modelle klären gleich viel Varianz auf}\) kann man mit dem anova()-Befehl:

anova(mod04_sleepstudy, mod03_sleepstudy)
## refitting model(s) with ML (instead of REML)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
mod03_sleepstudy 4 1802.079 1814.850 -897.0393 1794.079 NA NA NA
mod04_sleepstudy 6 1763.939 1783.097 -875.9697 1751.939 42.1393 2 0

Unterschied ist signifikant, also steigert der zusätzliche Random Slope die Modellpassung überzufällig.

Multi-level Modelle in MPlus

## Version:  1.1.0
## We work hard to write this free software. Please help us get credit by citing: 
## 
## Hallquist, M. N. & Wiley, J. F. (2018). MplusAutomation: An R Package for Facilitating Large-Scale Latent Variable Analyses in Mplus. Structural Equation Modeling, 25, 621-638. doi: 10.1080/10705511.2017.1402334.
## 
## -- see citation("MplusAutomation").
## When hashfilename = FALSE, writeData cannot be 'ifmissing', setting to 'always'
## The file(s)
##  'mod.dat' 
## currently exist(s) and will be overwritten

MPlus ist statistisch ultra-mächtig aber im vgl. zu sehr umständlich, da es keine Data Wrangling Funktionen hat. Daher muss man seine Multi-Level Daten zunächst mit einer anderen Software aufbereiten und als .dat-File speichern, bevor man dieses wieder in MPlus einliest, in einem .inp-Textfile die Syntax speichert und ausführt, aufdass die Ergebnisse dann in einem wiederum gesondeten .out-Textfile gespeichert werden. Da es im Mplus keine Objektorientierung gibt, muss man dann Modellvergleiche (z.B. Devianztest) händisch anhand der beiden .out-Files machen.

`

Die Syntax sieht in MPlus wie folgt aus:

TITLE:    Random Intercept Random Slope Modell (sleepstudy)
DATA:     FILE = "MPlus/mod.dat";
 
VARIABLE: NAMES ARE reac days subj;
          CLUSTER = subj;
          WITHIN = days;
ANALYSIS: TYPE = TWOLEVEL RANDOM;
MODEL:    %WITHIN%
          s | reac ON days;
          %BETWEEN%
          reac s;
          

JASP

In JASP/jamovi bekommt ihr das Random Intercept Random Slope Modell, wenn ihr den sleepstudy Datensatz als .sav-File einlest und dann die folgenden Einstellungen vornehmt:

Einstellungen für ein Random Intercept Random Slope Modell

Einstellungen für ein Random Intercept Random Slope Modell

oder gleich diese reproduzierbare JASP-Analyse öffnet.

Literatur

Hox, J. J., Moerbeek, M., & Schoot, R. van de. (2017). Multilevel analysis: Techniques and applications (Third edition). Routledge.

Gelman, A., & Hill, J. (2007). Data analysis using regression and multilevel/hierarchical models (Bd. 1). Cambridge University Press.

Snijders, T. A., & Bosker, R. J. (2012). Multilevel analysis: An introduction to basic and advanced multilevel modeling (2nd ed). Sage.